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Using many-body perturbation theory and coupled-cluster theory, we calculate the ground-state 
energy of *He and ^^O. We perform these calculations using a no-core G-matrix interaction derived 
from a realistic nucleon-nucleon potential. Our calculations employ up to two-particle-two-hole 
coupled-cluster amplitudes. 
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I. INTRODUCTION 



The coupled-cluster method originated in nuclear 
physics over forty years ago when Coester and Kummel 
proposed an exponential ansatz to describe correlations 
within a nucleus 0, Q- This ansatz has been well justi- 
fied for many-body problems using a formalism in which 
the cluster functions are constructed by cluster operators 
acting on a reference determinant Q . Early applications 
to finite nuclei were described in Ref. |3| . From that time 
to this, a systematic development and implementation of 
this interesting many-body theory in nuclear physics ap- 
plications has been only sporadic. The view from com- 
putational quantum chemistry is quite different. In fact, 
coupled-cluster methods applied to cornputational chem- 
istry enjoy tremendous success H, IE S S IS ^| over 
a broad class of chemistry problems related to chemical 
and molecular structure and chemical reactions. 

Many solid theoretical reasons exist that motivate a 
pursuit of coupled-cluster methods. First of all, the 
method is fully microscopic and is capable of system- 
atic and hierarchical improvements. Indeed, when one 
expands the cluster operator in coupled-cluster theory to 
all A particles in the system, one exactly produces the 
fully-correlated many-body wave function of the system. 
The only input that the method requires is the nucleon- 
nucleon interaction. The method may also be extended 
to higher-order interactions such as the three-nucleon in- 
teraction. Second, the method is size extensive which 
means that only linked diagrams appear in the computa- 
tion of the energy (the expectation value of the Hamilto- 
nian) and amplitude equations. As discussed in Ref. 
all shell model calculations that use particle-hole trunca- 
tion schemes actually suffer from the inclusion of uncon- 
nected diagrams in computations of the energy. Third, 
coupled-cluster theory is also size consistent which means 
that the energy of two non-interacting fragments com- 
puted separately is the same as that computed for both 
fragments simultaneously. In chemistry, where the study 
of reactions is quite important, this is a crucial property 
not available in the interacting shell model (named con- 
figuration interaction in chemistry). Fourth, while the 
theory is not variational, the energy behaves as a varia- 
tional quantity in most instances. Finally, from a com- 



putational point of view, the practical implementation of 
coupled cluster theory is amenable to parallel computing. 

Applications to nuclear problems resurfaced a few 
years ago in the works of Mihalia and Heisenberg 0, 0, 
Il3i . .14] . These efforts focused primarily on the structure 
of ^^O, and used a strategy of solution that is somewhat 
different from the approach we will take in this and sub- 
sequent articles. One major difference is that we will use 
a G-matrix to renormalize the two-body interactions be- 
fore we begin our coupled-cluster calculations. We will 
also take a somewhat different approach in our Hilbert 
space truncation. Also notable is the work of Moliner, 
Walet, and Bishop |(15|| who are pursuing nuclear prob- 
lems in translationally invariant coupled cluster methods 
in coordinate space. 

The computed energy using the coupled-cluster for- 
malism includes a very large class of many-body pertur- 
bation theory diagrams. In standard many-body pertur- 
bation theory, one typically sums all diagrams order by 
order. The coupled-cluster approach essentially iterates 
diagrams so that one may discuss it in terms of an infi- 
nite summation of particular classes of diagrams. Thus 
the theory is nonperturbative. In fact the coupled-cluster 
energy at the single and double excitation level contains 
contributions identical to those of second order and third 
order many-body perturbation theory, but lacks triple ex- 
citation contributions necessary to complete fourth-order 
many-body perturbation theory; see e.g., the review ar- 
ticle of Bartlett 5] . It has been shown that the quadru- 
ple excitation contributions may be factored exactly into 
products of double excitations, but no such factorization 
is possible for the corresponding triples. Therefore, the 
coupled-cluster energy lacks only triple excitation contri- 
butions to be complete through fourth order. 

In this paper, we wish to establish a line of research 
that we intend to pursue for calculating nuclear proper- 
ties using coupled-cluster techniques. This is therefore 
a first paper in a series that we will publish to both 
develop the method for nuclear physics and to demon- 
strate the power of the method for various applications. 
This first installment will be devoted to outlining our 
approach, investigating the physical motivations, estab- 
lishing numerical convergence tests, and presenting some 
initial calculations using the method. In Sec. ^ we will 
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describe our choice of reduced Hilbert space, construc- 
tion of an effective interaction for various model spaces 
and elimination of the spurious center of mass motion. 
The model spaces are defined in terms of various major 
harmonic oscillator shells and the effective interaction 
is defined in terms of the G-matrix, see e.g., Ref. [Tfij . 
Our calculations of the binding energy of nuclei like he- 
lium and oxygen entail therefore a dependence upon the 
number of harmonic oscillator states, the oscillator pa- 
rameter and the starting energy at which the G-matrix 
is computed. The convergence of the binding energy as 
function of the number of harmonic oscillator shells is a 
crucial test of the method discussed in this work. Since 
the coupled-cluster calculations are rather time consum- 
ing for systems like -'^^O, we present, as an introduction 
to the coupled-cluster method, intermediate results from 
perturbative many-body approaches in Sec. IIIII Sec. IIIII 
also serves the purpose of finding an eventual minimum 
for the energy as function of the oscillator parameter, 
with which we limit the number of coupled-cluster com- 
putations. It may also be of help in finding how many 
oscillator shells are needed in order to achieve a satisfac- 
tory convergence. In addition, this section sheds lights on 
the limitations of many-body perturbation theory com- 
pared with the coupled-cluster approach. 

We turn in Sec. IIVI to a description of the coupled- 
cluster equations. We discuss the numerical techniques 
we will employ to obtain solutions to the coupled-cluster 
equations and demonstrate several results for ^He and 
^^O in Sec.0 We conclude with a prospective for future 
directions of this research in Sec. IVII 



its nonperturbative low-energy limit at energy scales that 
are characteristic for low-energy nuclear physics. One 
promising way to circumvent this problem is to employ 
a derivation of the nuclear force based on chiral effective 
field theory jH 111- The authors of Ref. ^ under- 
took the task of generating an accurate nucleon-nucleon 
interaction based on chiral perturbation theory. They 
included one- and two-pion exchange contributions up 
to chiral order three. They also showed that a quan- 
titative fit of the nucleon-nucleon D-wave phase shifts 
requires contact terms representing short-range forces of 
order four. The number of free parameters used in this 
chiral interaction is 46, which is similar to the number of 
free parameters found in other two-nucleon forces. The 
phase-shift analysis shows excellent agreement between 
the chiral interaction and the scattering data. 

Two interactions were formulated in Ref. These 
two models, denoted Idaho-A and Idaho-B, differ in their 
prediction of the I?-state probabilities of the deuteron, 
while both interactions will give the same values for the 
■^S*!, '^Di and ei phase parameters up to 300 MeV in 
scattering energy. Idaho-A yields a D-state probability 
of 4.17%, while Idaho-B gives 4.94%. This also affects the 
triton binding energy, yielding 8.14 MeV and 8.02 MeV 
for Idaho-A and Idaho-B, respectively. A similar interac- 
tion which now goes to fourth order in chiral perturbation 
theory and includes charge dependence, has recently been 
presented by Entem and Machleidt, see Ref. Since 
this is a methodological paper, we limit the attention to 
one of these interactions, namely the Idaho-A model. Re- 
sults for other interactions, such as the Vis model of the 
Argonne group 'l8'|, will be presented in future work. 



II. EFFECTIVE INTERACTIONS FOR A 
TRUNCATED HILBERT SPACE 

The aim of this section is to present and partly justify 
the computation of an effective two-body Hamiltonian 
acting within a reduced Hilbert space. This two-body 
Hamiltonian will in turn serve as the starting point for 
the perturbative approach of Sec. IIIII and the coupled- 
cluster expansion discussed in Sec. IIVI 

Before we can compute such an effective two-body 
Hamiltonian, we need to define the nucleon-nucleon inter- 
action. Several types of modern nucleon-nucleon scatter- 
ing interactions have been developed during recent years. 
These interactions all fit nucleon-nucleon scattering data 
up to 300 MeV with excellent precision |il[ll[li|. They 
do give slightly differing results for the radius of the 
deuteron, the binding energy of the triton and also con- 
tain slight differences in the way they treat locality. 

Very recent work by Entem and Machleidt ^Oj pro- 
vides for the first time an interaction of quantitative ac- 
curacy that is based on effective field theory. One basic 
open question of nuclear theory involves understanding 
how the nucleon-nucleon interaction may be derived from 
quantum chromodymanics, the theory of strong interac- 
tions. Quantum chromodymanics has not been solved in 



A. Definition of the model space and tlie two-body 
effective interaction 

In order to derive an effective interaction suitable for 
coupled cluster calculations, we need to introduce vari- 
ous notations and definitions pertinent to the methods 
exposed. 

A common practice in nuclear many-body theory is 
to reduce the infinitely many degrees of freedom of the 
Hilbert space to those represented by a physically moti- 
vated subspace, the model space. In such truncations of 
the Hilbert space, the notions of a projection operator P 
onto the model space and its complement Q are intro- 
duced. The projection operators defining the model and 
excluded spaces are 

P = f]|$,)($,|, (1) 

and 

oo 
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with D being the dimension of the model space, and 
PQ ^ 0, ^ P, Q'^ = Q and P + Q = I. The two- 
body wave functions are normally eigenfunctions of 
an unperturbed Hamiltonian Hq- In this work we let 
only the kinetic energy enter the definition of Hq, i.e., 
Ho = t. Since we will employ a harmonic oscillator ba- 
sis, this means that we need to compute the expectation 
value of Hq as well. The unperturbed wave functions 
are not eigenfunctions of t. The full Hamiltonian is then 
H — t + VjvAT with V/vAT the nucleon-nucleon interaction. 
The eigenfunctions of the full two-body Hamiltonian are 
denoted by |^'q) and Ea, 

= S„ . (3) 

Rather than solving the full Schrodinger equation above, 
we define an effective Hamiltonian acting within the 
model space such that 

Pi/effP|*a) =^aP|*a) =Sa|$a) (4) 



where = P\'i'a) is the projection of the full wave 
function onto the model space, the model space wave 
function. Here H^s is an effective Hamiltonian acting 
solely within the chosen model space given by Hcs = 
PtP + Vcff, with the interaction 



where V^g\ V^^\ V^g\... are effective one-body, two- 
body, three-body interactions etc. For finite ^-body sys- 
tems, the sum terminates at i = A. As stated above, in 
this work we will limit the attention to two-body inter- 
actions. The next step could be to employ perturbative 
many-body techniques or the coupled cluster method. In 
perturbation theory, the effective interaction Hcff can be 
written out order by order in the interaction Vnn as 



PH,ffP = PtP + PVnnP + PVnn-VnnP + PVnn-Vnn-VnnP +■■■■ (6) 

e e e 

I 



In this expansion, e = uj — t, where ld is the so-called 
starting energy, defined as the unperturbed energy of the 
interacting particles. However, the nucleon-nucleon in- 
teractions all possess a hard core that makes them un- 
suitable for perturbative many-body approaches. The 
standard procedure is therefore to renormalize the short- 
range part of the interaction by introducing the so-called 
reaction matrix G 

G^Vnn + Vnn (7) 

Lo — QtQ 

The operator Q is normally different from the projec- 
tion operator defined in Eq. since the G-matrix by 
construction allows only specific two-body states to be 
defined by Q. Typically, the G-matrix is the sum over 
all ladder type of diagrams with intermediate particle- 
particle states only. This sum is meant to renormalize 
the repulsive short-range part of the interaction. The 
physical interpretation is that the particles must interact 
with each other an infinite number of times in order to 
produce a finite interaction. This interaction can in turn 
serve as an effective interaction acting in a reduced space. 

We illustrate the definition of the exclusion operator 
employed in this work in Fig. ^ Using a harmonic oscil- 
lator basis for the single-particle wave functions, a single- 
particle state is classified by the quantum numbers nlj. 
A two-particle state in an angular momentum coupling 
scheme is given by \(nalajani3lf}ji3)JTz) , where a and (3 
represent one of the orbitals Osi/2, 0p3/2, etc and J 

is the total two-particle angular momentum and Tz the 



corresponding isospin projection. 

The single-particle states labeled by niliji and 712^2^/2 
represent the last orbit of model space P. In this 
work niliji and 712/2^2 will mark the number of har- 
monic oscillator shells included in the definition of P. 
In the actual calculations presented below these range 
from four to eight major shells. For four major shells 
nihji = 2si/2 and n2?2j2 = 2si/2 while for eight 
major shells we get niliji = 3pi/2 and 712^2^2 = 
3pi/2 as the last single-particle orbits. In Fig. ^ the 
two-body state \{nalajan0lpji3)JTz) does not belong 
to the model space and is included in the computa- 
tion of the G-matrix. Similarly, \{nalajanjl^jj)JTz) 
and \{nslsjsnf3lf3jf3)JTz) also enter the definition of Q 
whereas \{nslsjs'n'yl^j-y)JTz) is not included in the com- 
putation of G. This means that correlations not defined 
in the G-matrix need to be computed by other non- 
perturbative resummations or many-body schemes. This 
is where the coupled-cluster scheme enters. 

Before we proceed we outline the computation of the 
G-matrix using the exclusion operator of Fig. ^ One 
can solve the equation for the G-matrix for finite nuclei 
by employing a formally exact technique for handling Q 
discussed in e.g., Ref. [ig- Using the matrix identity, for 
which P is the complement of Q such that P + Q = 1, 

QeQ e e Pe'^P e' ^ ' 
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nshjs 




FIG. 1: Definition of the exclusion operator used to compute 
the G-matrix. See text for further details. 



to rewrite Eq. Q as 

G = Gf + AG, , 
where Gi? is the free G-matrix defined as 



Of — Vnn + Vnn 



1 



UJ ~ t 



(9) 



(10) 



The term AG is a correction term defined entirely within 
the model space P and given by 

AG^^V^r^jPj^PjV^N. (11) 

Employing the definition for the free G-matrix of Eq. 
(|10|1 . one can rewrite the latter equation as 



AG= -Gf-P- 

e P(e-i 



~1 



—^P-Gf, (12) 



iGFe-i)P e 



with e = u! — t. We see then that the G-matrix is ex- 
pressed as the smTL of two terms; the first term is the free 
G-matrix with no corrections included, while the second 
term accounts for medium modifications due to the exclu- 
sion operator Q. The second term can easily be obtained 
by some simple matrix operations involving the model- 
space matrix P only. The above allows, for a given model 
space operator P, for a numerically exact computation of 
the G-matrix. 

The G-matrix defined by the exclusion operator Q of 
Fig. n represents now our effective two-body interaction, 
with a reduced Hilbert space defined by the model space 
P. This means that our Hamiltonian acts within the 
model space P and is given by 



The results from exact shell-model diagonalizations, per- 
turbative many-body calculations and coupled-cluster 
calculations will depend on the limits of the model space 
P and the chosen oscillator parameter. Furthermore, al- 
though the G-matrix has a weak dependence upon the 
starting energy oj, this dependence will affect the final 
results the calculations. The dependence of the results 
upon these parameters will be elucidated in Sees. IIIll and 

m 

We end this subsection with the setup of the exclusion 
operator used in the final perturbative many-body calcu- 
lations of Sec, mil and the coupled cluster calculations of 
Sec. El 

With the G-matrix model space P of Fig. we can 
now define an appropriate space for many-body pertur- 
bation theory and coupled-cluster calculations where cor- 
relations not included in the G-matrix are to be gener- 
ated. This model space is defined in Fig. [3 where the la- 
bels Up^qlp^qjp^q represent the same single-particle orbits 
as ni,2^i,2ji,2 in Fig. Hereafter we use the notation 
that p, q, r, s index all single-particle states, while i, j, fc, I 
refer to single-hole states 

The G-matrix does not refiect a specific nucleus and 
thereby single-particle orbits which define the uncorre- 
lated Slater determinant. For a nucleus like *He the 
Osi/2 orbit is fully occupied and defines thereby single- 
hole states. These are labeled by riiliji in Fig. [21 For 
^^O the corresponding hole states are represented by the 
orbits Osi/2, OP3/2 ^^d Opi/2- With this caveat we can 
then generate correlations not included in the G-matrix. 

The unperturbed Hamiltonian Hq is given by the ki- 
netic energy only. However, in order to define a suitable 
starting point for an effective interaction to be used in 
the coupled-cluster calculations we use the definition of 
particles and holes in Fig. Inland compute the single- hole 
energies through 



H^sH^t + Giio). 



(13) 



= {i\ + = + 1*^') ' (14) 

where F stands for the Fermi energy. We do not per- 
form a self-consistent Brueckner-Hartree-Fock calcula- 
tion however, as done by e.g.. Gad and Miither [23|, the 
reason being that this self-consistency is taken care of 
by the coupled-cluster procedure. The matrix elements 
are all antisymmetrized. Furthermore, for single-particle 
states above the Fermi energy we leave the single-particle 
energies unchanged. This procedure, which follows the 
Bethe-Brandow-Petschek theorem |25] , introduces an ar- 
tificial gap at the Fermi surface. Note also that the single- 
particle wave functions are not changed in Eq. H14|) . The 
main purpose of the above procedure is to yield a pre- 
scription for obtaining a starting energy independent ef- 
fective interaction for the coupled-cluster calculations. 
Using the single-particle energies from Eq. (|14|l we de- 
fine, following Ref. [l^l, an effective interaction for our 
coupled-cluster model spaces by 

1 



mv.s\ki) 
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FIG. 2: Definition of particle and hole states for coupled- 
cluster and perturbative many-body calculations. We use the 
notation that p, q, r, s index all single-particle states, while 
i,j,k,l refer to single-hole states The orbits represented by 
the quantum numbers riihii represent hole states whereas 
np,qlp,qjp,q represent the last particle orbits included in the 
G-matrix model space. The hole states define the Fermi en- 
ergy. 



+ {tj\G{oj = ek+ei)\kl)], (15) 
for two-body states with holes only and 

{ip\Vcff\kq) = ^[{ip\G{uj = e,+ep)\kq) 

+ {ip\G{uj = ek + eg)\kq)], (16) 

for two-body states with one particle and one hole. 
For two-body states with two single-particle states pq 
we use a fixed starting energy, typically in the range 
uj € [—80, —5] MeV. This introduces a starting energy 
dependence in our results. The reason for fixing the start- 
ing energies for two-particle states is due to the fact that 
we use kinetic energies only above the Fermi surface and 
our G-matrices are computed at negative starting ener- 
gies only. 

An obvious improvement to this procedure is to gen- 
erate a complex G-matrix which takes care of both pos- 
itive and negative energies, reflecting thereby the non- 
resonant continuum, bound two-body states and even- 
tual resonances and/or virtual states. We will however 
leave this interesting topic to a further study. The main 
focus of this paper is to establish a method for doing 
the coupled-cluster calculations without having to resort 
to employing the bare interaction, needing thereby many 
major harmonic oscillator shells. Our hope is that with a 
G-matrix, 7 — 8 major shells may suffice. Such a trunca- 
tion in the harmonic oscillator space is supported by the 



recent works of Barrett, Navratil, and Vary |26l 1271 128|. 
see also the recent calculations of Ref . |2^ . In these works 
no-core shell-model calculations have been mounted for 
light nuclei ranging from the triton to mass A = 12. 
Furthermore, in Refs. 26^ 27}, another approach for ob- 
taining a starting energy independent interaction is ob- 
tained through the similarity transformations of Lee and 
Suzuki 30, 31], yielding a fully hermitian effective inter- 
action. This should be contrasted to the more approxi- 
mative method presented in Eqs. H15|) and H16|l . 



B. Treatment of center-of-mass motion 

Momentum conservation requires that a many-body 
wave function must factorize as \l/(r) = 0(R)*(rrci) 
where R is the center-of-mass coordinate and Trei the 
relative coordinates. If we choose to expand our wave 
functions in the harmonic oscillator basis, then we are 
able to exactly separate the center-of-mass motion from 
the problem provided that we work in a model space that 
includes all nhuj excitations. In this paper, we perform 
calculations with a Q operator that allows for all possi- 
ble two-particle interactions within a given set of oscilla- 
tor shells. This means that we are nhto incomplete in a 
given calculation so that our method of separation of the 
center-of-mass motion becomes approximate. For exam- 
ple, for "^He in four major oscillator shells, we can excite 
all particles to n = 12?iti; excitations, but we can only ex- 
cite one particle to n — 3huj excitations. Thus, care must 
be taken when correcting for center-of-mass contamina- 
tion in our calculations. We have taken a variational ap- 
proach based on the the work of Whitehead et ai, |33 |. 
The idea is to add /3c.m.-f^c.m to the Hamiltonian, but with 
/3c. m. remaining fairly small. This minimizes the effects 
of the center-of-mass contamination on low-lying state 
properties, and partially pushes unwanted states out of 
the spectrum. If we were to use a large /?c.m., we would 
find spurious states entering into the calculated low-lying 
spectrum due to the incompleteness of our model space. 

We proceed as follows. The center-of-mass Hamilto- 
nian is 



1 



2MA 



-huj 



(17) 



where P = I]j=i,aP» ^ = iJ2i=i.A^i)/^- -^c.m. 
can be rewritten as a one-body harmonic potential, and 
a two-body term that depends on both the relative and 
center-of-mass coordinates of the two interacting parti- 
cles. The matrix elements for the two body terms may 
be found in Ref. [s^- Operationally, we add Hem to our 
Hamiltonian 



H' — H + (3c. m. Hem , 



(18) 



where we choose /3c. m so that the expectation value of 



Hr 



is zero l34l. This insures that our center-of-mass 



contamination within the many-body wave function is 
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minimized. We also find that this procedure yields rea- 
sonable spectra (in a space of four major oscillator shells) 
for ^He 1^3 . 

III. PERTURBATIVE MANY-BODY METHODS 

Our results will depend on the size of the model space 
and the chosen harmonic oscillator energy huj. This sec- 
tion serves therefore two aims central to the coupled- 
cluster approach of this work. 

• In the coupled-cluster calculations we search for an 
energy minimum as function of e.g., the oscilla- 
tor energy Tllu. Even for small systems like ^^O 
with five or six major shells included, these are 
major time-consuming calculations. Results from 
many-body perturbation theory, which to second 
or third order in the interaction G are fairly sim- 
ple, may therefore serve as a guide in order to limit 
the range of hu values used in the coupled-cluster 
calculations. With increasing dimensionality of the 
G-matrix model space, the results should become 
independent of the oscillator energy. 

• The energy will also depend on the size of the 
model space. The hope is that not too many shells 
are needed in order to achieve a converged energy. 
With the present approach, our coupled-cluster cal- 
culations of ^^O are limited to six major shells. Re- 
sults from many-body perturbation theory can thus 
serve as a guideline. 

The linked-diagram theorem 0, [s^ can be used to 
obtain a perturbative expansion for the energy in terms 
of the perturbation V{G) or V = H — Hq where Ho 
represents the unperturbed part of the Hamiltonian. The 
expression for the energy E reads 

oo 

E = Y,{^,\H[{cj- Ho)-'H] ' |*o)l , (19) 

fe=0 

where ^'o is the uncorrelated Slater determinant for the 
ground state, u is the corresponding unperturbed energy 
and the subscript L stands for linked diagrams only. In 



our calculations, we must replace the Hamiltonian H 
with the effective one defined in Eq. H13|) and employ 
the definition of particle and hole states of Fig. |21 In 
Fig. |31 we show all antisymmetrized Goldstone diagrams 
through third order in perturbation theory (we omit the 
first order diagram). All closed circles stand for a sum- 
mation over hole states. In this section we let the G- 
matrix define the interaction vertex. Since we do not 
use a self-consistently determined single-particle Hamil- 
tonian, we need to account for diagrams with so-called 
Hartree-Fock insertions as well. Examples of the latter 
are shown in Fig. O see also the work of Kassis [s^ for 
a detailed discussion of the various diagrams. There is 
an additional problem with our many-body perturbation 




13 14 15 16 



FIG. 3; Antisymmetrized Goldstone diagrams through third 
order in perturbation theory included in the evaluation of the 
binding energy. The dashed lines represents the interaction, 
in our case the G-matrix. Particle and hole states are rep- 
resented by upward and downward arrows, respectively. The 
first order diagram is omitted. All closed circles stand for a 
summation over hole states. 



theory calculations. Consider the expression AE^ for di- 
agram 4 in an angular momentum coupled basis (with J 
being the total two-body angular momentum and the 
corresponding isospin projection) of Fig. |21 



= ^ {'2J+i)mJT.\G{uj = e,+e,)\ipq)JT,)— ^ (20) 

o — Si + £j — £p ~ Sq 

ij < F 
pqrs > F 
J 

X {{pq)JT,\ Giuj = £, + £j) \{rs)JT,) — — {{rs)JT,\G{u = + e,) . 

> ^ j £r ^ s 

I 



The G-matrices depend on the starting energy w, defined cillator basis and with kinetic energies only these start- 
in this case by the hole energies. With a harmonic os- 
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ing energies will be positive, whereas our G-matrix is 
defined for negative energies only with lo £ [—140,-5] 
MeV. To remove this problem we employ the single-hole 
energies computed according to Eq. H14|l. This leads to 
an artificial gap at the Fermi surface and demonstrates 
one of the problems with many-body perturbation the- 
ory. A self-consistent approach is needed both for holes 
and particles, and the G-matrix needs to be computed 
for both positive and negative energies. This is however 
beyond the scope of this work, where our emphasis is on 
the coupled-cluster approach. The center-of-mass correc- 
tions discussed in the previous section are not included 
in our perturbative calculations. This section, as will be 
seen below, serves the main aim of justifying our model 
space used in the coupled-cluster computation. In ad- 
dition, the energy minimum as function of the oscillator 
energy limits the range of Tllo values used in the coupled- 
cluster calculations. 

Needless to say, many-body perturbation theory has 
severe limitations. It is very difficult to go beyond third 
order in perturbation theory without a self-consistent 
single-particle potential and beyond fourth order with 
a self-consistent potential. The coupled-cluster method 
offers on the other hand a systematic way to generate all 
many-body correlations in a given model space. 



A. Many-body perturbation theory results for 
helium and oxygen 

We present here results from third-order in perturba- 
tion theory for the binding energies of "^He and "'^^O as 
functions of the size of the model space and the chosen 
oscillator energy hu). These results are shown in Figs. 01 
and[Slfor ''He and ^^O, respectively. 





FIG. 4: Binding energy E from third-order perturbation the- 
ory for ''He as function of the number of major harmonic 
oscillator shells N and the oscillator energy Tilo. For N = 8 
we have the optimal value of E — —20.830 MeV at hu) = 13.3 
MeV. The experimental value is E = -28 MeV. 



hm (MeV) 

FIG. 5: Binding energy per particle E/A from third-order 
perturbation theory for as function of the number of ma- 
jor harmonic oscillator shells N and the oscillator energy Tiu. 
For iV = 8 we have the optimal value of E/A = -7.120 MeV 
at huj = 13.6 MeV. The experimental value is E/A = -7.98 
MeV. 



There are several features to be noted. First of all, 
both figures show that the results seem to stabilize be- 
tween seven and eight major shells. For ''He all possible 
excitations within these shells are allowed in the compu- 
tation of the diagrams, to be contrasted to the work of 
Kassis 37] and other traditional approaches 16] where 
one typically considers only 2 — Ahw excitations. For '^^O 
we keep the same numbers of maximum allowed excita- 
tions for both the Osi/2 shell and the OP3/2OP1/2 shells. 
As an example, for eight major shells we could have lAtko 
2p — 2h excitations from the Osi/2 shell whereas from the 
Op3/20pi/2 shell we can at most have 12hoj 2p — 2h ex- 
citations. The latter fixes the total number of allowed 
excitations with eight major shells for '^^O. The fact that 
the energies seem to converge at this level of truncation is 
a welcome feature which can be exploited in the coupled- 
cluster calculations. These calculations, see below, are 
much more challenging from a computational point of 
view since we in principle generate a much larger class of 
diagrams. In this work we are limited to computations 
with at most seven major shells in helium and six major 
shells in oxygen in our coupled-cluster calculations. This 
means that hopefully the trend seen in Figs. 0] and fal- 
lows us to limit our coupled-cluster calculations to six 
or seven major shells. Secondly, although the minimum 
shifts a little as function of the oscillator energy as we in- 
crease the oscillator space, we notice that as the number 
of major shells is increased, the dependence of the bind- 
ing energy upon the oscillator parameter weakens. A 
similar feature is seen in the coupled-cluster calculations 
below. For '^^O the minimum for seven shells takes place 
at E/A = -7.155 MeV for tuv = 12.9 MeV and for eight 
shells we have E/A = -7.120 MeV at huj = 13.6 MeV. 
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For six shells we obtain EjA = -7.416 MeV at huj = 12.6 
MeV. The curvature for larger values of fiLO decreases 
with increasing number of shells N . At fiuo — 18 MeV 
we have for ^^O and TV = 8 that d{E/A)/dijj 0.217, 
for iV = 7 we obtain d{E/A)/duj = 0.277 and N ^ Q we 
have d{E / A) / dui — 0.346. The corresponding derivatives 
for ^He at hcj = 18 MeV are for iV = 8 d{E/A)/duj = 
0.119, for TV = 7 d{E/A)/dw = 0.158 and for = 6 
d{E/A)/du} = 0.209, indicating a smoother dependence 
upon TiLO with increasing TV. The reader should also note 
that in the limit huj — > we have i? — > 0. Finally, we 
observe that we are not able to reproduce the experimen- 
tal binding energies. We will come back to this point in 
Sec. where a comparison with coupled cluster theory 
is also made. 

I 



where Eq = {'^q \ Hcff{uj) \ 5'o) and the Fock-matrix 
element is given by 

fp, = {p\t\q)+Y,{P^\\ri) . (22) 

i 

To simplify notation, we will use {pq \ \ rs) ~ {pq \ G{uj) \ 
rs). We employ bracket notation {} to indicate normal 
ordering with respect to the reference state. As stated 
previously, we use the notation that p, q, r, s refers to all 
single-particle states and i,j,k,l index all sums below 
the fermi surface. In addition we let a, b, c, d index all 
sums above the fermi surface. The total number of single- 
particle states in the model space is Ns = Np + Nh where 
Np refers to the number of particle states, and TV;i is the 
number of hole states. Hereafter we simply use H for 
ifeff(w)- 

Formally, the coupled-cluster method begins by pos- 
tulating that the correlated many-body wave function is 
given by 

= exp (T) I *o) , (23) 

where we define the the correlation operator as 

T = Ti + + Ta + • ■ ■ + Ta . (24) 

The correlation operators are defined in terms of n- 
particle n-hole (np-nh) excitation amplitudes as 

Ti ^ *>t^^ . (25) 

i<6f,a>ef 

T2 = itj^tata^a, , (26) 

i,j<e f;ab>ef 



IV. COUPLED CLUSTERS IN SINGLE AND 
DOUBLE EXCITATIONS 

In this section, we will discuss a formal derivation 
of the coupled-cluster equations. While this discus- 
sion is standard in quantum chemistry (see, for exam- 
ple Ref. 6]), the nuclear physics community may not be 
familiar with this formulation of coupled-cluster theory. 
We will therefore introduce notations and equations that 
we will continue to use throughout our discussion (in this 
and following papers) of this technique and its extensions. 

We must first bring the Hamiltonian into normal or- 
dered form with respect to the reference state | "^o)- The 
Hamiltonian then becomes 



(21) 

I 

and higher order terms for T^, to Ta- Coupled-cluster 
theory may thus be hierarchically improved upon by in- 
creasing the number of Ti operators one computes. We 
will call the theory in which only Ti and T2 operators are 
present, CCSD, or coupled-clusters at the single and dou- 
ble excitation level. CCSDT means that T3 is retained in 
the correlation operator, while CCSDTQ refers to keep- 
ing both T3 and T4 correlation operators. In this un- 
coupled representation, the correlation amplitudes must 
obey the fermion-symmetry relations which for the T2 
correlation operators yield tfj = —tjf = —t\'- = t''^. We 
will use the short-hand notation ti and t2 to represent 
the array of all Ip-lh and 2p-2h operators. 

We compute the expectation value of the energy from 

E={^>o\ exp i-T) H exp (T) | *o) • (27) 

Because the energy is computed using projective, asym- 
metric techniques, an important question concerns the 
physical reality of the coupled-cluster energy. Quantum 
mechanics requires that physical observables should be 
expectation values of Hermitian operators. The coupled- 
cluster energy expression contains the non-Herniitian op- 
erator [exp(— T)ff exp(T)]. However, if T is not trun- 
cated, the similarity-transformed operator exhibits an 
energy-eigenvalue spectrum that is identical to the origi- 
nal Hermitian operator, thus justifying its formal use. 
From a practical point of view, the coupled-cluster energy 
tends to follow the expectation value result (if the theory 
is reformulated as a variational theory), even when T is 
truncated. 

The correlation energy is quite easy to calculate and is 
given by 



Hcs{cj) = t + G'(w) = ^/pq{a],a,} + i ^(pg II rs){aj;ajasar} 

pq pqrs 

+ {^Q I H,s{uj) I ^-o) =Hn + Eo, 
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ia aibj aibj 



For two-body Hamiltonians, this equation is general and 
is not restricted to the CCSD approximation since higher- 
order cluster operators such as and T4 cannot produce 
fully contracted terms with the Hamiltonian and there- 
fore contribute zero to the energy. Higher-order opera- 
tors can contribute to the energy indirectly through the 
equations used to determine these amplitudes. The three 
terms in Eq. H28|) are usually referred to as the Ti , T2 , and 
Ti contributions to the correlation energy. 

The equations for amplitudes are found by left projec- 
tion of excited Slater determinants so that 

= (vI/nexp(-r)i/^(r)|vI/o), (29) 
= {^f^ I exp (-T) Hn (T) | M-o) . (30) 

The Baker-Hausdorf relation may be used to rewrite the 
similarity transformation as 

exp (-T) Hn (T) = Hn + [Hn, Ti] + [iJ^, T2] 

-fi[[i7Ar,Ti],Ti] + i[[i7jv,T2],T2] 

+ [[Hn,Ti],T2] + -- - . (31) 



This equation is non-linear in the ti amplitudes, and lin- 
ear in the t2 amplitudes. 



I 

The expansion terminates exactly at quadruply nested 
commutators when the Hamiltonian contains at most 
two-body terms, and at six nested commutators when 
three-body terms are present. We stress that this termi- 
nation is exact, thus allowing for a derivation of exact 
expressions for the ti and t2 amplitudes. To derive these 
equations is straightforward but tedious work. For ex- 
ample, one of the amplitude terms is given by 

pqrs kc Id mnef 

(*o I {afajabaajia^a'^asar} 

{a^ak}{a+ak}{a^ai}{ata'j:anarn} \ *o) , (32) 

for which Wick's theorem may be used to calculate the 
expectation matrix element. Here Gat is the normal or- 
dered G-matrix operator, and the subscript c indicates 
only connected diagrams enter into the computation of 
the expectation value. 

The ti amplitude equations are given by 



(33) 



The t2 amplitude equations are given by 



I 

= fa^+Yl fact^ - E ^-^^fc + E(^« 1 1 + E /^-^^fc + ^ E(^« 1 1 

c h kc kc kcd 

- ^ E(^^ II c*)*^? - E /''-^^'^^ - E(^^ II c*)*^^" + E<^« II ^^)^^^' 

klc kc klc kcd 

- ^(fcz II cd)tititi + ^(fcni ^d)titt -\Y.^^i II cd)tttt - ^ E<^^ II "^^^ ■ 

klcd klcd klcd klcd 

I 



J 



= (a6 II zj) + ^ (ifc,t- - /,,t^) - ^ (/fe,t^ - Uafk) 

c k 

+ \ Y.^kl II rj)tl\ + \ Y.'^ab II cd)<^^ + P{ij)P{ab)Y,{kb \\ cj)C 

kl cd kc 

+ P{ij)Y.{ab II cM - P{ab)Y,{kb II ij)tl 

c k 

+ \p{i3)P{ab)Y,{kl II cd)t1^tf^ + J E(^^ II - \n(^b)Y.^kl II cd)t1^tfi - \p{ij) ^(fc^l 

klcd klcd klcd klcd 
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+ ^Pi'^b) Y^ikl II + ip(zj) Y^iab II cd)e,t'^ - P{tj)P{ab) Y,{kh \\ ic)tlf^ 

kl cd kc 

+ P{ab) J2 hct^ + E f'^'^^ 

kc kc 

+ P{i3)P{ah) Y,{kl\\ ic)t1t% + ipfe-) Y.^kl\\ ck)i!itt - ip(a5) ^(fc6 || cd)tlt<-^ 

klc klc kcd 

- \p{ij)P{ah)Y,{kh II cdntltl + \p(i3)P{ah)Y,{kl II cj)t^t^tf 

- P(i3)Y.(kl II cd)i=,tfi?,^ - P{ah) Y^lkl II cd)t^fetrt^^ + ip(^j-) ^(fc; II cd)t,=t^^t^l' 

/c^crf A;^crf klcd 

+ ip(a6) ^ (fc; 1 1 + P{i3)P{oh) Y^(kl\\ cd)t1t\tlj + -P{i3)P{ah) ^ {kl \ \ cd)t1tltjt 



(34) 



klcd 



klcd 



klcd 



The permutation operator P yields 

P{i])f{i]) - /(*J) - /(ii) (35) 

The equations of the ^2 ampHtudes are nonhnear in both 
ti and t2 terms. While these equations appear quite 
lengthy, they are solvable through iterative techniques 
that we will discuss below. We note that the ampli- 
tude equations include terms that allow for 4p-4h ex- 
citations. Indeed, while we speak of doubles in terms of 
amplitudes, the class of diagrams involved in the theory 
include fourth-order terms. This is a very important dif- 
ference and distinction between the shell model with up 
to 2p-2h excitations and CCSD. Furthermore, when the 
energy is computed in CCSD, all terms are linked and 
connected. 

In order to calculate expectation values of operators 
we use the Hellmann-Feynman theorem [3^ which states 
that if we perturb our Hamiltonian such that H' ^ H + 
Xn where A is a small quantity and O is the operator of 
interest, then the energy changes only by a small amount 
from its original value of E{X ~ 0). As a function of A, 
the energy becomes E' = E{X = 0) -|- XdE/dX, and the 
expectation value of the operator is given by 

V. CCSD CALCULATIONS FOR HELIUM AND 
OXYGEN 

A. Iteration of the equations 

Several computational challenges arise when we imple- 
ment the CCSD equation solver. One problem involves 



memory requirements for the G-matrix, which in its un- 
coupled form is a 4-index tensor and therefore requires a 
large amount of storage. In order to maintain fast com- 
putation, we do not employ storage compression of the 
G-matrix at this stage. Thus, for example, an iV = 7 cal- 
culation requires 100 GBytes of storage for the G-matrix 
elements. Present-day parallel computing systems use 
distributed memory architectures that allow for the stor- 
age of such large sets of data. Our implementation in 
distributing the G-matrix is to store the third and fourth 
indices across processors in sub-matrix blocks. In doing 
this, we are able to take advantage of the large available 
memory for storing the matrix elements. 

We use an iterative method to generate solutions to 
the ti and t2 amplitude equations. In this paper, we are 
concerned with closed shell nuclei. In this case, a slight 
rearrangement of Eqs. and (|^ gives a second order 
perturbative initial solution for the amplitudes. For the 
ti amplitude, we rearrange the first few terms by pulling 
out the diagonal in the first two sums and defining a 
Ip-lh and 2p-2h energy denominators as 



Dt = U - faa , (37) 
P'tj = fii + fjj - faa - fbb ■ (38) 



The first terms in the ti and t2 amplitude equations then 
become 
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Dtt^ = fa^++ 5^(1 - S,a)fact', " ^(1 " S.k)^^ + ' ' ' , 
c k 

D^t^ = (a6|Uj)+P(a6)^(l-5fc,)A,t- + P(y)^(l-4,)/.,t?f + 



(39) 
(40) 



While these are exactly the same equations as given in 
Eqs. and we may use them to begin an iterative 
solution by initially setting all amplitudes on the left- 
hand side to zero. We then obtain for initial amplitudes 



,a _ Jai 



j.ab 



fa 

ab 1 1 ij) 



(41) 
(42) 



We now need to compute the various terms in Eq. H34() 
to obtain the new amplitudes. Since our G- matrix ele- 
ments are distributed across processors, we will have par- 
tial sums on certain indices. For clarity, we do not show 
this complication in the following. We demonstrate the 
numerical procedure by considering one of the terms in 
the two-particle-two- hole amplitude equation, Eq. H34|) : 



fiab,ij) 



E 

kl.cd 



{kl \\ cd)t'rftt1 



(43) 



We perform the following mapping to matrices M, TV, 
and O, P, and Q: 



{kl II cd) 

cd 

f{ab,ij) 



t. 



(44) 



where we map the indices a — (a, 6), /3 = {k,l), 7 = 
(c, d), and S — Our computation then becomes 

two matrix-matrix multiplications: 



Pa,7 = E ^-^"./S^; 
Qa..5 — ^ ^ ^a.-yO. 



(45) 



followed by a mapping of Qa,5 to f{ab,ij). All terms 
in the amplitude equations can be formulated in this 
way and three important consequences follow. First, the 
computational work scales from one computation of or- 
der 0{N^N^) to two computations of 0{N^Nl). Sec- 
ond, by formulating the algorithm in this manner, we 
take advantage of high ly o ptimized Basic Linear Algebra 
Subprograms (BLAS) [3g|. This represents a significant 
reduction in effort when contrasted with the naive equa- 
tions. We see that a calculation ^^O should require six- 
teen times more computational time when compared to 



N Ns ^He 


G (in millions) 


16 Q 


4 80 1,792 


0.93 


24,960 


5 140 4,000 


7.23 


77,880 


6 224 7,976 


40.4 


176,240 


7 336 14,112 


178.5 





TABLE L Various memory scalings for the He and O sys- 
tems as a function of the number of major oscillator shells, 
N. Listed are the number of uncoupled single-particle states, 
the two-particle-two-hole amplitudes in *He, the number of 
non-zero G-matrix elements, and the number of nonzero two- 
particle-two-hole amplitudes in ^^O. 



a calculation of ^He in the same model space. Due to the 
longer loop structure in ^^0 the time required to com- 
plete a single oxygen run is only a factor of ten larger than 
the He case. Third, because we have sub-blocked the in- 
teraction amplitudes, we can spread the work across nu- 
merous processors. In doing so, we must perform a global 
reduction operation (a global sum) in order to generate 
the new amplitudes for the next iteration. We show in 
Table some details of the computational sizes of the 
m-scheme problems that we have undertaken in this pa- 
per. The number of unknowns for which we are solving is 
approximately the number of two-particle-two-hole am- 
plitudes. Therefore, in our largest calculation for ^^O we 
are actually solving approximately 176,300 equations per 
iteration. We also list the number of nonzero G-matrix 
elements in Tabled 

Eq. 142|) represents the terms that one uses to com- 
pute the energy in second-order perturbation theory of 
the M0ller-Plesset type We show in Table HH the 

energy obtained from the 0'^ £^(0) iteration and the final 
iteration E{F) as a function of increasing oscillator levels 
in the ^^O system. Note that the difference between the 
converged CCSD energies and the initial O"* order ener- 
gies increases as the basis space increases. The converged 
summation of the CCSD equations yields approximately 
10 MeV (or 0.6 MeV per particle) in extra binding. These 
findings are corroborated by those from many-body per- 
turbation theory from Sec. IIIII It is therefore worth 
comparing these results with those from second-order 
and third-order many-body perturbation theory as well. 
These are labeled £^mbpt ^-nd i?MBPT the same table. 
The reader should notice that the zeroth iterations of 
the coupled-cluster schemes already includes corrections 
to the one-body amplitudes ti. However, the energy de- 
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4 14 -135.117 -132.063 -129.920 -140.473 

5 14 -124.786 -124.844 -121.520 -127.790 

6 14 -121.356 -121.481 -118.233 -119.733 



TABLE II: Comparisons of the order energy E{0) and 
the converged CCSD results E{F) for /3c.m = 0.0 in ^^O as a 
function of increasing model model space. The results are also 
compared with many-body perturbation theory to second and 
third order, S^^px and ^mbpti respectively. All energies are 
in MeV. 



nominators used in the computation of the second-order 
diagrams of Fig. O (diagrams 2 and 3) have hole states 
determined by Eq. (|14|l . The agreement with the zcroth 
order iteration and second-order perturbation theory is 
very good, especially for five and six major shells, as 
can be seen from Tabled However, for third-order per- 
turbation theory one clearly sees fairly large differences 
compared with the coupled-cluster results. Typically, the 
relation between first- and second-order in perturbation 
theory for -'^^O is given by a factor of ^ 5 — 6. For e.g., 
TV = 5 and = 14 MeV, we have --329.124 MeV from 
first order and —47.719 from second order. To third or- 
der we obtain a repulsive contribution of 3.324 MeV, to 
be contrasted with the almost 3 MeV of attraction given 
by higher-order terms in the coupled-cluster expansion. 
This indicates that many-body perturbation theory to 
third order is most likely not a converged result. An 
interesting feature to be noted from many-body pertur- 
bation theory calculations is that higher terms loose their 
importance as the size of system is increased. For ^He 
the relation between first-order and second-order pertur- 
bation theory is given by a factor of ^ 3 — 4, depending 
on the value of hcj. Calculations for ^"Ca not reported 
here indicate a relation of ~ 7 — 9 between first-order and 
second-order perturbation theory. This is somewhat ex- 
pected since the G-matrix is smaller for larger systems, 
although the energy denominators becomes smaller. 

We show typical convergence curves for our coupled- 
cluster calculations in Fig. El for ^^O for the N = 4,5,6 
oscillator shell model spaces. In this figure, E{I) is the 
energy at the iteration, while E{F) is the final en- 
ergy of the system. We obtain convergence within 50 
iterations in most cases. The curves also exhibit an ex- 
ponentially convergent behavior. 

Since the CCSD amplitude equations are nonlinear, we 
investigated the results obtained by using different initial 
reference states. The naive choice for | ^o) is to fill the 
lowest oscillator states for a given system. For example, 
we define | '^q) as the filled Os state for ^He and the filled 
Os-Op states for ^^O. This choice is fine for closed shell 
nuclei since none of the energy denominators discussed 
above becomes zero, but for other nuclei this would be a 
problem. As an alternative procedure, we start with the 
Hartree-Fock ground-state Slater determinant for a given 



> 



' 1 ' 1 




1 ' 


■ N ^ ■, 

^ 

\ • ■ 




HF basis, N=4 

-- OSC basis, N=4 








■-■ OSC basis, N=5 








■-■ OSC basis, N=6 




. - - 

\^ \ ^. 


















r 






1 , 1 




' s 

^. 1 





I I I I I I I L l_ 

10 20 30 40 



CCSD Iteration 

FIG. 6: The convergence of the ground-state energy as a func- 
tion of the CCSD iterations for ^''O. 



system. We solve the Hartree-Fock equations in the oscil- 
lator basis in order to obtain transformation matrices D 
that take us from the oscillator to the Hartree-Fock basis. 
We then transform the Hamiltonian to the Hartree-Fock 
basis using the relation aj = J2a Daicl^ where ai and a\ 
annihilate and create particles in the oscillator basis and 
Cc and cj, annihilate and create particles in the Hartree- 
Fock basis. Note that D = 1. While a complete di- 
agonalization of H in either basis would yield the same 
results, the CCSD amplitude equations are not invariant 
under this transformation since states below and above 
the fermi surface will be mixed. Furthermore, at the 
Hartree-Fock level, rotational symmetry of the Hamil- 
tonian is broken, although correlations including those 
at the CCSD level will restore much of this symmetry. 
Although we do not discuss in this paper open-shell sys- 
tems, the Hartree-Fock solution offers a clean way obtain 
a reference Slater determinant for those cases. 

We show in Fig.Elthe convergence of the CCSD equa- 
tions in both the oscillator and Hartree-Fock basis for 
= 4 major shells. In both cases, convergence at the 
10~^ level is reached before 40 iterations. A more qual- 
itative assessment of how the choice of reference state 
affects the final results is seen in Table IIIII In the os- 
cillator basis, the Ip — \h amplitudes carry a significant 
fraction of the correlation energy, while in the Hartree- 
Fock basis these terms do not contribute to the energy. 
The difference between the Hartree-Fock and oscillator 
total energies is 0.6%. We note that if one completely 
ignores the Ip ~ \h amplitudes in the Hartree-Fock ba- 
sis, the energy becomes -138.38 MeV, a difference of less 
than 0.07%. In the oscillator basis, a much larger error 
of 5.6% is obtained if one ignores the Ip — Ih ampli- 
tudes. Furthermore, the CCD equations (ignoring the 
Ip — Ih amplitudes) sometimes exhibit numerical insta- 
bilities similar to ignoring the Hartree-Fock insertions in 
many-body perturbation theory discussed above. 
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Term Oscillator (MeV) Hart 


ree-Fock (MeV) 


Eo 


-109.45 


-121.977 


Ti 


-9.669 


7 X 10"" 




-1.757 


-0.3 X 10"^ 


T2 


-18.439 


-16.495 




-29.865 


-16.498 




-139.31 


-138.47 



TABLE III: Comparisons of CCSD results in ^"^O when us- 
ing naively filled oscillator reference state, or when using 
the Hartree-Fock reference state. These calculations were per- 
formed at hu = 14 Mev and /3c. m. = 1.0. 
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FIG. 7: Dependence of the ground-state energy of *He on Tuj) 
as a function of increasing model space. 



B. ''He and ^"O ground states 

We now return to a discussion of *He and ^^O by pro- 
viding a description of their ground-state energies using 
the CCSD formalism. As in the many-body perturbation 
theory section of this paper, we wish to demonstrate how 
the coupled-cluster theory converges as a function of in- 
creasing model space. We are currently able to perform 
this study for up to seven major oscillator shells in helium 
and up to six shells in oxygen. While these studies do not 
address the starting-energy dependence of the G matrix, 
they do indicate the convergence of the calculations as 
a function of model space, and they indicate the soften- 
ing of the hio dependence as one moves to larger spaces. 
The fiuj dependence of the ground state energy of ^He 
as a function of the model space is shown in Fig. for 
the TV = 4, 5,6, 7 major oscillator shells with /3c. m. = 0. 
These curves generally exhibit a parabolic character. 

We applied the center-of-mass correction described 
above to both the He and O nuclei. We demonstrate 
how this procedure behaves when one solves the CCSD 
equations in Fig. |S1 for ''He in the TV = 4 and model 




FIG. 8: Top panel: the total energy of ''He in four major 
oscillator shells plotted as a function of the center-of-mass 
parameter /3c. m.. Bottom panel: the center of mass energy as 
a function of /3c. m. . The energy is at a variational minimum in 
this 4 oscillator shell space at luo = 10 MeV. The variational 
minimum for as a function of the CCSD iterations. 



space. We calculated the expectation value of of iJc.m. 
with Eq. H36|) . These plots indicate that the ccntcr-of- 
mass behavior is less severe when we perform calculations 
at the variational minimum in ?iaj, which in this case is 
at Tiio = 10 MeV. We observe that these corrections are 
rather small for both He and O, amounting to a correc- 
tion of less than 1% in both cases. This correction only 
applies to the ground state. 

We performed a full diagonalization of ^He in four ma- 
jor oscillator shells at Tiu) = 10 MeV and /3c. m. = 0.5 and 
obtained an energy of —23.4837 MeV for the ground state 
as compared to our CCSD result of -21.98 MeV, a 6.8% 
difference. We anticipate that triples corrections to the 
energy |4ll| will recover much of the remaining difference 
between the complete result and the CCSD result. The 
CCSD results should be approximately of similar quality 
in larger model spaces so that we can expect an approx- 
imate complete diagonalization energy of —22.4 MeV in 
seven oscillator shells as compared to —21.0 MeV for the 
CCSD result. The ground-state energy using Idaho- A 
was quo ted as -27.4 MeV by Navratil and Ormand in 
Ref. YM- The major difference between our projected 
full diagonalization result and -27.4 MeV is most likely 
due to the differing treatments of the G-matrix. 

In this initial study we performed calculations of the 
^^O ground state for up to six major oscillator shells as 
a function of fiLO. Fig. |51 indicates the level of conver- 
gence of the energy per particle for TV = 4, 5, 6 shells. 
The experimental value resides at 7.98 MeV per parti- 
cle. While this calculation is not completely converged, 
it does show signs of convergence as the space is in- 
creased. By six oscillator shells, the fiu) dependence be- 
comes rather minimal. We find a ground-state binding 
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FIG. 9: Dependence of the ground-state energy of ^^O on Tuu 
as a function of increasing model space. 



energy of 7.6 MeV per particle in oxygen using the Idaho- 
A potential. Since the Coulomb interaction should give 
approximately 1 MeV/A of repulsion, and is not included 
in this calculation, we actually obtain approximately 6.6 
MeV of nuclear binding in the 6 major shell calculation 
which is somewhat above the experimental value. We 
note that the entire procedure (G-matrix plus CCSD) 
tends to approach from below converged solutions. 



VI. CONCLUSIONS AND PERSPECTIVES 

Our goal in this paper has been to describe a nonper- 
turbative method of solution to the many-body problem 
that sums classes of diagrams built upon low order (up 
to fourth order in this paper) many-body perturbation 
theory diagrams. We have shown how to calculate nu- 
clear ground states using coupled-cluster methods. In 
this paper, we concentrated on CCSD equations. We 
used a G-matrix as our two-body interaction and we ex- 
panded in the spherical harmonic-oscillator basis. We 
reiterate that CCSD is a nonperturbative approach to 
the many-body problem: we sum classes of diagrams to 
infinite order in order to calculate various nuclear prop- 
erties. The coupled-cluster method discussed here clearly 
demonstrates the need of summing many-body correla- 
tions to infinite order. This is seen when comparing the 
coupled-cluster results to many-body perturbation the- 
ory. Furthermore, the latter is hard to extend beyond 
third order without a self-consistent single-particle po- 
tential and beyond fourth order with a self-consistent 
potential. 

Before closing this paper, we would like to discuss sev- 
eral steps that we will take during the course of this re- 
search. 

One improvement upon the method will be to include 
the calculation of triples excitations (called CCSDT or 



approximations to it). We indicated that the three- 
particle-three-hole diagrams likely give repulsion and are 
important for the description of ground-state properties. 
Since the CCSD equations do not include the 3p-3h dia- 
grams completely, and since we have seen that these dia- 
grams are important, we will eventually need to include 
the triples amplitudes into our equations. Various meth- 
ods that include triples diagrams have been investigated 
by quantum chemists and we will investigate which of 
these methods are appropriate for the nuclear problem. 

The effort to perform a complete solution to the quan- 
tum many-body problem grows exponentially as one adds 
particles. Combined with the difficulty of methods we 
suffer in nuclear physics from interactions that are not 
completely determined. Thus, our methods and tech- 
niques for solution typically develop in lock step with 
our understanding of the nuclear Hamiltonian. While 
several sets of two-nucleon of interactions that fit nu- 
cleon scattering perfectly have been developed over the 
last 10 years, their many-body characteristics (in partic- 
ular, their ability to obtain nuclear ground-state masses) 
indicate that they are insufficient. Three- nucleon inter- 
actions become necessary even to fit the triton and "^He. 
To date, no derivations of CCSD (or CCSDT) equations 
exist that incorporate a three-body interaction. We will 
pursue this effort in future research. 

CCSD and its extensions can be used to ob- 
tain excited state information by diagonalizing H — 
exp (— T) H exp (T) in the space of all singly- and doubly- 
excited determinants where the amplitudes are obtained 
directly from the converged CCSD amplitudes. This will 
be an important step in the development of the coupled- 
cluster method for nuclear science. 



Finally, our results do depend on the starting energy 
of the effective interaction. For ^^O and four oscilla- 
tor shells and hu! = 14 MeV, our /3c. m. = 0.0 result 
is -140.47 MeV with the starting energy of -80 MeV, 
while we obtain —143.53 MeV with a starting energy of 
—60 MeV. The dependence is therefore weak, but still 
present. The dependence is more crucial in helium since 
the binding energy is much lower, and our G-matrix is 
not defined for positive starting energies. There are two 
possible ways to overcome this problem. One is to use 
the similarity transformation of Lee and Suzuki, follow- 
ing closely the no-core approach of Barrett, Navratil and 
co-workers 26^ ^J. This yields a hermitian and start- 
ing energy independent interaction for a large space. Al- 
ternatively, one can compute self-consistently the single- 
particle energies using a G-matrix defined for both pos- 
itive and negative starting energies. Thereafter, a start- 
ing energy independent interaction can be obtained using 
e.g., the prescription of Eq. (|15|l . but for both holes and 
particles. These approaches will be investigated in future 
works. 
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